Intraspecific variation in a predator changes intertidal community through effects on a foundation species

Abstract Intraspecific variation is an important form of biodiversity that can alter community and ecosystem properties. Recent work demonstrates the community effects of intraspecific variation in predators via altering prey communities and in foundation species via shaping habitat attributes. However, tests of the community effects of intraspecific trait variation in predators acting on foundation species are lacking despite the fact that consumption of foundation species can have strong community effects by shaping habitat structure. Here, we tested the hypothesis that intraspecific foraging differences among populations of mussel‐drilling dogwhelk predators (Nucella) differentially alter intertidal communities through effects on foundational mussels. We conducted a 9‐month field experiment where we exposed intertidal mussel bed communities to predation from three Nucella populations that exhibit differences in size‐selectivity and consumption time for mussel prey. At the end of the experiment, we measured mussel bed structure, species diversity, and community composition. While exposure to Nucella originating from different populations did not significantly alter overall community diversity, we found that differences in Nucella mussel selectivity significantly altered foundational mussel bed structure, which in turn altered the biomass of shore crabs and periwinkle snails. Our study extends the emerging paradigm of the ecological importance of intraspecific variation to include the effects of intraspecific variation on predators of foundation species.

Rocky intertidal communities are often shaped by foundational mussels. In the California Current Ecosystem, the abundance and size structure of intertidal California mussels (Mytilus californianus) shape the community of coexisting species. The introduction of predators on mussels, either sea stars (Pisaster) or dogwhelks (Nucella), reduces the abundance of mussels, not only creating opportunities for competing species but also decreasing habitat for colonization and refuge for commensal species like small gastropods, worms, and crustaceans (Lafferty & Suchanek, 2016;Menge et al., 2021;Navarrete, 1996). While prior work has focused on the community effects of mussel predator presence or absence (Menge et al., 1994;Paine, 1966), there is a lack of understanding about the community effects of mussel predator intraspecific trait variation.
Nucella ostrina and Nucella emarginata are sister species of intertidal dogwhelks commonly found on the west coast of North America ( Figure 1). These species live in mid-intertidal California mussel beds and feed primarily on mussels and barnacles by drilling (Clelland & Saleuddin, 2000;West, 1986). In central California, their distributions overlap, and the two species are largely indistinguishable genetically, morphologically, and ecologically (Marko, 1998, 2005, Marko et al., 2014, including in a recent neutral genetic analysis of the populations used in the present study (Contolini, Reid, & Palkovacs, 2020). Therefore, we consider trait differences among these three populations equivalent to intraspecific trait variation and hereafter refer to the species collectively as "Nucella." Nucella have crawl-away larvae that develop and hatch in the same intertidal area as their parents. This life history limits gene flow and can lead to local adaptation in foraging traits (Dawson et al., 2014;Marko, 1998;Sanford & Worth, 2009). Indeed, Nucella exhibits population-level differences in size selectivity and consumption rate of mussels, which are traits that could alter the size structure of mussel beds (Contolini, Kroeker, & Palkovacs, 2020;Contolini, Reid, & Palkovacs, 2020).
Here we used a 9-month field experiment to test the prediction that population-level variation in Nucella consumption rate and size selectivity of mussel prey indirectly alters intertidal community composition by changing the physical structure of the mussel bed on which other organisms depend for inhabitable space. We predict that Nucella populations will differentially alter mussel density and size, and these changes will cascade to affect species diversity and community composition of organisms living within and atop the mussel bed matrix. Our study tests the role of intraspecific variation in predator foraging in communities through effects on foundation species.  Contolini, Reid, & Palkovacs, 2020). Nucella from Lompoc on average drill larger mussels (mean ± SE 54.18 ± 8.73 mm mussel length compared with 41.88 ± 1.35 and 42.13 ± 2.87 mm for Hopkins and Soberanes, respectively) and can consume a California mussel several days faster than Nucella from both Hopkins and Soberanes (mean ± SE 25.52 ± 2.92 d compared with 31.10 ± 2.68 and 27.52 ± 3.27 d; Contolini, Reid, & Palkovacs, 2020, Contolini, Kroeker, & Palkovacs, 2020. Nucella were held in filtered, flowing F I G U R E 1 Nucella dogwhelk (a) resting and (b) feeding with proboscis inserted into a California mussel.

| Mussel bed predation experiment
To test the effects of population-level variation in Nucella consumption rate and size selectivity of mussels, we performed a common garden experiment at Terrace Point in Santa Cruz, CA, USA (36.94°, −122.06°; Figure 2). Terrace Point is an exposed, gently sloping, southfacing site with a substrate of consolidated mudstone and multiple rocky benches separated by narrow sandy beaches. In the mid zone, the California mussels are the dominant primary space-holding species, forming a monoculture covering tens of square meters in a single layer on the horizontal surfaces. The mussel bed creates habitat for mobile animals such as gastropods, crabs, and worms; sedentary and sessile animals such as encrusting worms and barnacles; and algae including scouring pad algae (Endocladia muricata), Turkish towel (Mastocarpus spp.), and red turf algae (Calacanthus sp.). Terrace Point on average has seawater temperatures greater than that of Soberanes and Lompoc and lower than Hopkins, and of the four sites, it has a higher mean and lower standard deviation pH (Table S1).
In the experiment, we outplanted Nucella from all three populations to an experimental array of cages and let them feed within their assigned cages for 9 months. The cage array had a total area of 3.4 m 2 on an existing bench of continuous, level California mussel bed in the mid zone ( Figure S1). Mean Nucella shell length per cage initially was 24.22 ± 0.64 mm (mean ± SD; N = 24) and did not differ significantly among populations (ANOVA F 2,21 = 0.006, p = .99; Figure S2a). Mean Nucella mass per cage initially was 2.34 ± 0.24 g and also did not differ significantly among populations (Kruskal-Wallis test, X 2 21 = 23, p = .3; Table S2; Figure S2b). Any Nucella that died or were lost were replaced with one of a similar size from the same source population held in tanks in the nearby marine lab. We created the array by clearing all biological material 10 cm around thirty-two 20 × 20 cm natural mussel bed plots such that each was surrounded by a border of bare rock. We used existing mussel bed so we could test the responses of fully established mussel bed communities; for this reason, it was not possible to measure community diversity and composition within each cage at the start of the experiment to obtain a baseline because doing so required destructive sampling. We installed cages made from 0.4 cm Vexar mesh and secured them to the rock using stainless steel lag screws, washers, and marine epoxy (Z-Spar splash zone compound). Cages were 20 × 20 × 8 cm with removable lids secured with cable ties. We arranged the cages in rows parallel to the shore and assigned them to one of four treatments: five adult Nucella from one of the three populations or no Nucella (control), each replicated eight times and in a randomized block design ( Figure S3). Blocks (N = 8) ensured treatments were approximately equally exposed to edge and tidal effects and allowed us to account for this natural variation in statistical models. Nucella density was within the range of densities observed at the source sites (<1 to 10 m −2 ; Multi-Agency Rocky Intertidal Network (MARINe), 2021, pers. obs.). We marked Nucella with bright nail polish and uniquely numbered all individuals with bee tags (Bee Works, Canada) attached to the shell with cyanoacrylate glue. To control for differences in Nucella size, we measured Nucella starting and ending shell length from the shell apex to the tip of the anterior canal using digital calipers. We measured Nucella starting and ending mass by drying snails thoroughly with paper towels, squeezing excess water out from behind the operculum, and using an analytical scale (Mettler Toledo AG104). We initiated the experiment on 17 October 2017 by adding tagged and premeasured Nucella to assigned cages. We opened all cages biweekly to photograph the surface (Canon PowerShot D30), record drilled mussels, collect dislodged mussels, remove invading non-experimental Nucella, and replace dead or lost Nucella with one of a similar size from the same population (23 from Soberanes, 5 each from Hopkins and Lompoc). We monitored temperature (water and air) at the site and inside two cages every 15 min using HOBO temperature loggers (Onset Computer Corporation; Table S3) and compared this with F I G U R E 2 Map of study sites and experiment location with mean ± SD temperature and pH. Santa Cruz temperature is from a mid-zone intertidal logger for the duration of the experiment (17 Oct 2017 to 1 Aug 2018 [Multi-Agency Rocky Intertidal Network (MARINe) et al., 2017,2018]). Hopkins, Soberanes, and Santa Cruz pH data are from intertidal pH sensors from 15 July to 22 September. Lompoc data are from an offshore sensor (Purissima) from the same date range in 2011 (Rivest et al., 2016). More details of the pH and temperature regimes can be found in Table S1 and Chan et al., (2017).
The average mid zone temperature (air and water combined) at the experiment site was 13.31 ± 1.68 SD°C (Table S3). Inside cages, the mean combined air and water temperature was on average 0.52°C higher than the temperature of the mid-zone mussel bed outside of cages. Nine months later, on 29 July-1 August 2018, we measured the Nucella in the cages and removed and froze all mussels and their associated communities, including all organisms within and on top of the mussel bed matrix, at −20°C for further analyses.

| Mussel bed structure and community composition
We characterized the direct effects of Nucella on the mussel beds at the end of the experiment by measuring the sizes of drilled and remaining mussels using digital calipers. To characterize the ecological community within the mussel bed matrix, we identified all organisms within the mussel beds to the lowest possible taxon (https:// seanet.stanf ord.edu/; Light & Smith, 2007). We sorted all organisms, cleaned them of sand and large debris by hand or by rinsing them with fresh water, dried them in a 56°C oven (Chicago Surgical & Electrical Co. Imperial II, Thelco Precision Model 2, or Quincy Lab Model 40GC) for 7 d, and measured the dry mass of each taxon per cage using an analytical scale (Mettler Toledo AG104). Mass of all shelled organisms included shells. As we were unable to obtain accurate measurements of algal biomass due to it fracturing into small pieces when handled, we instead calculated the final percent cover of algae on top of the mussel beds, where most algae were located, using the bed surface photographs and ImageJ image analysis software v. 1.51 s (Abràmoff et al., 2004).

| Statistical analyses
We used analysis of variance (ANOVA) or, when data transformation did not normalize non-normal error distributions, Kruskal-Wallis tests to test for significant differences in Nucella sizes, number of drilled and remaining mussels, mean size of drilled and remaining mussels, community diversity (Shannon-Wiener diversity), and biomass of individual taxa among Nucella treatments. We used ANOVA models with type II sums of squares for number and size of drilled and remaining mussels, with main effects Nucella population and block to account for variation among blocks. We transformed variables when necessary to meet assumptions of normality and homogeneity of variances. We also performed permutational multivariate analysis of variance (PERMANOVA) to test for significant effects of Nucella population on community composition using Bray-Curtis dissimilarity community matrices, accounting for block as a main effect, and using 999 permutations. We used similarity percentages (SIMPER) with Bray-Curtis dissimilarities to find which species most discriminated between treatment communities (package vegan; Oksanen et al., 2018).
We tested for indirect effects of Nucella on diversity through changes in the mussel bed structure first using ANOVA to test for significant differences in the number and size of drilled and remaining mussels among treatments. We then used linear regressions to test for significant relationships between the number and size of drilled or remaining mussels and Shannon-Wiener diversity or biomass of influential taxa, as identified using SIMPER. Finally, we used   Table S4; Figure S2b).

| Direct effects of Nucella on mussel bed structure
Over the course of the experiment, Nucella from each population drilled over 150 mussels (N = 622 total) and visibly deteriorated mussel beds (Table S5). There were no significant differences in the number of mussels drilled among Nucella populations (ANOVA, block F 7,14 = 0.50, p = .82; population F 2,14 = 1.05, p = .30). The control cages had a low number of drilled mussels (N = 30) from small, transient native dogwhelk invaders that were immediately removed once found (mean ± SD 7.00 ± 4.10 invaders per 2-week sample period, and 133 total over the course of 9 months). The control treatment had significantly fewer drilled mussels than Nucella treatments (ANOVA, block F 7,21 = 0.47, p = .84; treatment F 3,21 = 14.63, p < .001; Tukey HSD p < .01 for populations paired with control).

| DISCUSS ION
The community effects of intraspecific trait variation are well known for predators shaping prey communities and for bottom-up effects of foundation species, but not for predators of foundation species.
In this study, we tested the effects of mussel predator trait variation on the mussel bed community. We used a 9-month field experiment to test the hypothesis that population-level variation in Nucella consumption rate and size selectivity of mussel prey indirectly alters intertidal community composition by changing the physical structure of the mussel bed on which other organisms depend. Our results supported this hypothesis. First, our results showed Nucella predation significantly affected mussel bed structure, as treatments with Nucella had significantly more and larger drilled mussels than the control treatment. Next, by drilling different sizes of mussels, Nucella populations differentially altered the structure of the mussel bed and the biomass of certain species living in the mussel bed matrix. These results provide the first evidence of the ecological effects of intraspecific variation in predators of foundation species.
Foraging adaptations could explain the variation in mean drilled mussel size among Nucella source populations we observed in our experiment. The process of attacking and consuming mussels is time-intensive and can take well over a week, during which time Nucella is exposed to abiotic and biotic stressors which can act as selective forces (Burrows & Hughes, 1989;Contolini, Kroeker, & Palkovacs, 2020). Nucella populations (including N. emarginata, N. ostrina, and Nucella canaliculata) throughout the California Current System show pronounced differences in foraging behaviors that have been linked with abiotic and biotic drivers (Contolini, Reid, & Palkovacs, 2020;Sanford & Worth, 2009).
For example, Soberanes and Lompoc naturally experience cooler temperatures than Hopkins (Figure 2; Chan et al., 2017). Under controlled conditions, Nucella from Soberanes and Lompoc consume mussels faster than those from Hopkins (Contolini, Kroeker, & Palkovacs, 2020), which is consistent with a pattern of countergradient variation, where elevated feeding rates maintain growth in cooler conditions (Yamahira & Conover, 2002). The mean temperature in our experiment was 2-3°C warmer than at Soberanes and Lompoc, but less than half a degree warmer than at Hopkins, TA B L E 1 Summary of ANOVA model of mean drilled mussel length versus Nucella population. and counter-gradient feeding adaptations could explain why the cooler source populations attacked larger, more caloric prey and had increased growth.
Our model showed evidence that Nucella affected shore crabs indirectly through habitat modification. Nucella population itself may not be the strongest direct driver of shore crab biomass because Nucella and shore crabs are not known to have strong direct ecological interactions. Mean drilled mussel length, however, may be a stronger direct driver of shore crab biomass because shore crabs use cracks and crevices for living spaces (Bovbjerg, 1960;Kanter, 1977). In our experiment, cages with larger drilled mussels had larger cracks and crevices and higher shore crab biomass, and we observed crabs inhabiting large empty mussel shells ( Figure S4). However, mean shore crab biomass was highest in the control treatment where mean shore crab size was smallest, implying that Nucella presence had a negative effect on juvenile shore crabs but a positive effect on larger adult shore crabs that increased with drilled mussel size.

Littorina spp. (periwinkle snails) were also indirectly affected by
Nucella through their interaction with shore crabs. Periwinkle snails are small-bodied (maximum shell length about 25 mm), and prey for shore crabs, so it follows that their biomass would be less in treatments where shore crab biomass was higher (Boulding et al., 2020).
Periwinkle biomass was not simply a product of mussel bed size structure because the relationship between periwinkle biomass and mean drilled mussel length was not significant when shore crabs were also included in the model. This is strong evidence that habitat modification by Nucella caused significant changes to the trophic interaction between shore crabs and periwinkles that depended on the population-specific pattern of Nucella predation on mussels.
Other studies on population-specific consumption patterns of marine consumers also confirm the importance of trait variation but do not often test for broader community effects. The presence or absence of salt marsh predators had differential effects on the diversity and biomass of two plant species by amplifying trait variation in herbivores (Hughes et al., 2015) and differences in Pisaster body condition were related to mussel cover (Menge et al., 2021 formal analysis (lead); funding acquisition (equal); investigation (equal); methodology (lead); project administration (equal); resources (equal); validation (equal); visualization (lead); writing -original draft (lead); writing -review and editing (equal). Eric P. Palkovacs: Conceptualization (equal); formal analysis (supporting); funding acquisition (equal); investigation (equal); methodology (supporting); project administration (equal); resources (equal); supervision (lead); validation (equal); visualization (equal); writing -original draft (supporting); writing -review and editing (equal).

ACK N OWLED G M ENTS
This study is part of G.C.'s Ph.D. dissertation research at the University of California Santa Cruz. We thank P. Raimondi, P.
Marko, E. Sanford, and E. Grosholz for helpful comments that improved the manuscript. S. Traverso, X. Clare, M. Ansaldi, E. White, and K. Stuhldreher provided invaluable help in the field and lab.
The Partnership for the Interdisciplinary Studies of Coastal Oceans (PISCO) provided temperature and pH data. Climate.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflicts of interest.

O PE N R E S E A RCH BA D G E S
This article has earned an Open Data badge for making publicly available the digitally-shareable data necessary to reproduce the reported results. The data is available at [10.7291/D1XM4Q].

F I G U R E 4
Path diagram showing results of the structural equation model for Nucella population, mean drilled mussel length, mean Pachygrapsus crassipes biomass, and mean Littorina spp. biomass in the mussel bed. Standardized regression coefficients are shown above paths and significances are shown below. Red arrows indicate negative relationships. Nonsignificant paths are gray.

F I G U R E 5
Schematic showing effects of intraspecific variation in Nucella predation on Pachygrapsus crassipes (shore crabs) and Littorina spp. (periwinkle snails). Size selectivity of mussel prey differed for each Nucella population. This differential consumption led to differences in shore crab biomass, in turn altering periwinkle snail biomass.